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ABSTRACT 

A  class  of  raethods  for  the  numerical  Integration  of 
the  time  dependent,  one  velocity  transport  equation  in 
spherical  geometry  is  presented.   All  of  the  methods  pro- 
ceed by  differencing  first  in  the  angular  variable.   The 
resulting  differential-difference  equations  form  a  hyper- 
bolic system  which,  by  proper  centering,  is  in  a  well- 
known  normal  form.   These  equations  are  then  transformed 
into  two  systems  of  integral  equations  from  which  the 
final  difference  equations  are  obtained. 

An  iterative  procedure  for  solving  one  of  the  systems 
of  difference  equations  is  shown  to  converge  under  stated 
conditions  on  the  mesh  spacing.   The  difference  scheme 
for  the  hyperbolic  system  is  shown  to  be  unconditionally 
stable  and  convergent.   Assuming  that  the  solution  of 
this  hyperbolic  system  converges,  with  refinement  of  the 
angular  mesh,  to  the  solution  of  the  transport  equation, 
a  sufficient  condition  is  obtained  for  convergence  of  the 
numerical  solution  to  the  transport  solution. 

The  methods  and  analysis  presented  are  easily  extended 
to  many-velocity  and  plane  geometry  problems.  Steady  state 
problems  may  be  treated  by  introducing  obvious  simplifications, 
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1»  Introduction 

Many  problems  concerned  with  radiation  or  neutron 
transport  require  the  solution  of  the  linearized  BoLtzmann 
equation  (for  the  density,  !D(  t,r ,  iv  ,v) ,   of  photons  or  neu- 
trons )  in  the  general  form 


(1.0) 


j  ^  4^  +  ;^  •  "V  +  <5(t,r,v)  I  !D  =  S(a);t,^,w  jv)   . 


Here   r   is  a  position  vector,   v   is  the  particle  speed, 
tJ  is  a  unit  vector  in  the  direction  of  particle  velocity, 
(5(t,r,v)   is  the  total  cross  section  at   (t,r)   for  par- 
tides  with  speed  v,   and  S(  !D;t,r,  <»/ ;v)   represents  the 
birth  of  particles  at   (t,r)   with  speed  v  in  direction 
w  due  to  scattering,  emission,  fission  or  any  combination 
of  these  processes.   Analytic  solutions  of   (1.0)   with 
appropriate  boundary  and  initial  conditions  and  specific 
forms  of   S  are  known  only  for  extremely  special  cases. 
Thus,  numerical  procedures  for  approximating  solutions  of 
these  problems  have  been  devised  and  tested  with  encouraging 
results  [1,2,5].   Most  of  these  procedures  at  present  lack  the 
theoretical  justification  necessary  for  accepted  numerical 
methods,  that  is,  proof  of  the  convergence  and  stability 
of  the  difference  equations.   Furthermore,  those  problems 
In  which  S  contains  a  scattering  or  fission  integral 
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usually  Involve  an  iterative  solution  of  the  difference 
equations,  the  convergence  of  which  is  rarely  demonstrated. 
In  the  present  report  we  formulate  a  class  of  numer- 
ical procedures,  for  a  plane  or  spherically  symmetric  form 
of   (1.0),   of  which  one  is  essentially  the   S  -  method 
devised  by  B.  Carlson  [l].   Some  variations  of  the  basic 
method  are  indicated  which  are  Intended  to  better  approxi- 
mate different  physical  situations  (i.e.  strong  absorbtion, 
rapid  spatial  variation  of  absorbtion,  etc.).   In  the 
latter  part  of  the  report  we  attempt  a  theoretical  justi- 
fication of  a  particular  one  of  the  class  of  difference 
equations  presented.   These  results  are  partially  success- 
ful  in  that  sufficient  conditions  for  the  convergence  of 
the  iterative  solution  are  obtained  and  "stability  and  con- 
vergence" of  the  space-time  integration  for  a  fixed  angular 
mesh  is  demonstrated.   The  latter  result  is  an  unconditional 
convergence,  as  /\t  — >  0  and  /\r  — >  0,   to  the  solution 
of  a  system  of  hyperbolic  differential  equations.   'ATiat  re- 
mains to  be  proven  is  the  convergence,  as  the  angular  spac- 
ing goes  to  zero  (perhaps  not  independently  of  ^t   and 
/\r) .   of  the  solution  of  the  hyperbolic  system  to  the 
solution  of  the  Eoltzmann  equation. 

2.   Formulation 

We  consider  the  H.inear  Eoltzmann  equation  in  spherically 
symmetric  geometry  for  one  velocity  group 

-  5  - 


Here   !D  Is  a  particle  density,  C    a  total  cross  section 
for  the  loss  of  particles,   jx  is  the  cosine  of  the  angle 
between  the  radius  vector  and  velocity  vector,   v   is  the 
particle  speed  (also  a  measiire  of  its  energy)  and   S   is 
a  source  term.   We  take,  for  the  present, 

1 
(2.1a)  S|2);t,r,ii]  =  g(t,r,n)  +  v^^(r)|  (  !D( t,r,tt)dia  + 


-1 


+  <^^^v)\   (  K(ii,,tx')E(t,r,ii')dii'   , 


where 


(2.1b)  I  (  K(n.,ti')dti'  =  1 

-1 


The  inclusion  of  the  inhomogeneous  term,   g(t,r,ti.),   allows 
the  treatment  of  many-velocity  problems  within  the  present 
formulation. 

The  boundary  conditions,  for  simplicity,  assume  no 
particles  entering  the  sphere   0  <  r  <  R  from  r  >  R  : 

(2.2a)  a)(t,R,tA)  =  0  ,   11  <  0   . 
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For  continuity  in  angle  at  the  origin  we  require' 
(2.2b)  a}(t,0,ti)  =  iD(t,0,-ii)   . 

The  initial  particle  density  is  specified  by 


(2.3)  !D(0,r,ti)  =  iD°(r,ti) 


where  the  given  function  !D   satisfies   (2.2). 

For  infinite  plane  geometries  the  above  formulation 
becomes  valid  by  eliminating  the  —  term  in   (2.0),   and 
conditions  {2m2)      imply  a  symmetric  configuration  in  the 
Interval   -  R  <  x  <  R.   More  general  plane  cases  are  easily 
treated  by  replacing   (2.2b)   by  (D(t,0,ix)  =  0,  }x  >  0, 
Steady  state  problems  in  either  geometry  are  included  by 
letting  V  -->  00  in   (2.0)   and  dropping  the  initial 
condition   (2,3). 

3»  Approxiriiating  Hyperbolic  System 

Let  the  [x- Interval,  -1  <  IJ-  <  1,  be  divided  into 
J  Intervals  of  equal  length  ,  ^[i  =  2/ J,  by  the  mesh 
points 

(3.0)  tij  =  -1  +  i^[i   ,      0  <  j  <  J   . 

In  fact  for  such  spherically  symmetric  problems   ffi(t,0,^) 
must  be  Independent  of  [i,      but  we  use  (2.2b)  for  simplicity. 

'The  restriction  of  equal  mesh  spacing  could  easily  be  elim- 
inated and  is  Included  only  for  simplicity  of  formulation  and 
notation.   A  better  choice  of  the  \i ,   Is,  perhaps,  to  take  them 
at  the  zeros  of   Pj^,(n,),  «' 
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We  seek  functions   (/.(t,r)   which  approximate  the  solution 
(D(t,r,ti..)   on  this  mesh.   Equations  for  the  determination 
of  the   y  .   are  obtained  by  replacing  the   (x-derivatlve 

J 

in   (2.0)   by  the  difference  quotient   (  cj?  .-  rG.  -,  )/Ap.. 
centered  at  [i .   ^  /^  =   [i .-  /\[i./2,      replacing  the  other 
terms  by  approximations  centered  at  the  same  point  and 
multiplying  by  two.   This  procedure  yields  equations  of 
the  form 

The  quantities  a.,  a.,  b.   and  S.(t,r)   are  determined 
such  that  the  approximations 

(3.2b)      (^4|)  -ill  (cp  .-cp.  ,)   , 

(3.2c)      s(IDJt,r,^l]  *  i(S.+  S.  -) 

are  properly  centered  and  have  a  minimum  truncation  error 
when  the  q?  .(t,r)   are  replaced  by  the  \!l{t,r,[i.).      Some 
possible  choices,  all  having  truncations  of  the  same  order 
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''■'''         '^3^'  =.     ,  *  i<^j -^  ^  ^- -^^  '   ' 


and 


in     /Xtii      are   listed  for   the     a.      and     a.      In  table    I  and 

J  J 

for  the  b,   in  table  II.   For  simplicity  In  formulation 
and  analysis  we  define  the   S.  by  means  of  the  trapezoidal 
rule  as  follows: 


{3.3a) 


S.(t,r)  3  gj(t,r)  +  y   a^j.(r)  (p^(t,r)   , 


where 


(3.3b) 


g^(t,r)  =  g(t,r,tij) 


a^,(r)  =  [v<Tp(r)  +  cTgCr  )K(m.  ..ix^)]  ^  e^    ° 


^J=l 


e^  =  2,  i  5^  0,J   . 


J-^o 

A 

B 

c 

D 

"3 

¥^i^^3-l^ 

¥^^3*^i-l^ 

^j 

Jl    ^"^t^J^t^j-l^ 

^J 

I'^'j^^j-l' 

I'^'J^^^J-l' 

^j-i 

m+1    ^^^j'-^^'^j-l^ 

TABLE  I 


i£o 


B 


-^[l-(^ 


t^^t^.1l^,--l-*-^^tl)] 


2  2 
Ati^     2    ^^ 


TABLE  II 

The  entries  In  columns  B  of  the  tables  are  those  used 
in  [1]  and  are  determineQ  there  by  assuming  linear  dependence 


-  9  - 


on  \x     within  ti-intervals  and  Integrating  across  each 

interval.   The  arbitrariness  shovm  in  the  above  tables 

is  thus  lost.  We  prefer  the  coefficients  in  columns  A 

which  have  the  property  that  a.=  a.,   and  the  b.  are 

"exact"  in  the  sense  that  the  approximations   (3»2b)   are 
centered  at   n-.,  z^*   With  this  choice  equations   (3»1) 

can  be  written  as 

1  <  j  <  J     . 

These  are   J  equations  for  the   J+1   functions  ^.. 
Another  equation  is  obtained,  as  in  [l],  by  setting 
ti  =  u  =  -1   in   (2.0)   and  making  approximations  as 
above  to  get 

Equations   (3.i|.),  the  differential-difference  equa- 
tions which  approximate  the  transport  equation   (2.0), 
form  a  hyperbolic  system  of  first  order  linear  partial 
differential  equations.   Furthermore,  by  virtue  of  the 
choice  of  coefficients  for  which  a.  =  a.,   these  equa- 

J        J 

tions  are  in  a  well-known  normal  form  [3]  in  which  all 
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differentiated  variables  in  each  equation  are  differen- 
tiated in  the  same  direction.   The  importance  of  this  form 
will  become  apparent  in  the  "derivation"  of  the  difference 
equations  in  section  I|.  and  in  the  convergence  proof  of 
section  6. 

The  boundary  conditions,   (2,2),   are  approximated 
as  follows: 

(3.5a)  fj(t,R)  =  0   ,   J  <  J/2   ; 

(3.5b)  fj^**^^  =  fj-j(**0)   • 

The  initial  conditions   (2.3)   are  replaced  by 

(3.6)  <?j(0,r)  =  iD°(r,tij)   . 

It  can  be  shown  that  for  sufficiently  smooth  3)  (r,tJ,.) 
the  above  initial  and  boundary  conditions  are  sufficient 
to  determine  a  unique  solution  of  the  hyperbolic  system 
(3.i^.)(see  [ij.]). 

The  transport  problem  has  now  been  replaced  by  an 
approximating  initial-boundary  value  problem  for  a  hyper- 
bolic system  of  equations.   As  an  aid  in  the  derivation 
of  the  final  difference  equations  it  is  useful  to  consida:' 
some  other  forms  of   (3.l4.)« 
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Since  these  equations  are  in  the  normal  form  we  may 
replace  them  by  an  equivalent  system  of  integral  equations, 
as  is  usual  in  the  theory  of  hyperbolic  systems,  by  inte- 
grating along  the  characteristics.   The  characteristics 
are  straight  lines  in  the   (t,r)   plane,  defined  by 


, ,  „  V    dt  _  1     dr  _  li    T^  _  /  2  .  T  /  2   . 


where  £,.     represents  arc-length  along  the  j-th  charac- 
teristic.  Equations   (3»i+a)   can  be  written  as 


Equation   (3«4-b)   may  also  be  written  in  the  above  form  by 
defining 

(3.9)  a^  =  -1,  b^  =  0,  (f_^  =  S_^  =  0  . 


Integrating   (3.8)   from   P.   to  Q.   along  the  j-th  charac- 
teristic yields  the  system  of  integral  equations 

Q 

(3.10)  (fj-fj.i)^  =  (fj.fj.,)^.pSj.Sj.,.(^bjA)fj- 


0      J       J"-^  p 


J 


J 
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Another  useful  system  of  integral  equations  is  ob- 
tained by  writing   (3»8)   as 

(3.11)  D^  af]"fj^^3-l'  *   (^-''jAXfj^fj.l)  =  Sj.  Sj.^  * 

Then  assuming  the  r.h.s.  knovm  we  solve  each  of   (3«11) 
as  an  ordinary  differential  equation  along  the  appropriate 

characteristic  to  obtain         q 

n        b  ,   <^^1 
'^   (c3+^j/r)  -pi- 

(3.12)  (f^.-^^.i)  =  (?j-fj.i)   e   J-  ^   ^ 

^JV^.I^^^M-^^        ^^ 


U..   Finite  Difference  Equations 

We  shall  "derive"  the  final  difference  equations  by 
placing  a  mesh; 

r^  =   k/\r,  Ar  s  ^^K  ,   0  <  k  <  K  ; 

(ij..O) 

t^  =  n^t ,  At  S  /N  ,   0  <  n  <  N   ; 

-  13  - 


In  the  rectangle   (0<r<R,  0<t<T)   and  numerically 

integrating  the  system  (3«i|-)   along  characteristics  on 

this  mesh.   This  is  essentially  the  procedure  of 

Carlson  [l].   Thus  in  (3-10)  we  identify  all  the   Q. 

with  an  arbitrary  point  Q,  =    (t  ^,  ,r,  )   of  the  mesh   (I4..O) 

The   P.   are  taken  to  be  the  first  intersection  of  a  mesh 
J 

line  with  the  j-th  back  characteristic  through  Q.   As 
the  characteristics,  defined  in   (3.7),  are  straight 
lines  with  slopes 

(i|.l)  Xj.  =  Vvaj  ,   0  <  j  <  J   , 

the  points   P.   are  conveniently  classified  into  four 
groups  as  follows  (see  Pig.  1.):  Let  the  ratio  of  the 
mesh  slope,  ^-^   //\r,   to  the  characteristic  slope  be 
designated  by 


(^•2)  6^  =^^  ,   0<  j  <  J 


Then  a  point   P.   shall  belong  to  one  of  the  classes  A  , 
A_,  B_|.  or  B__  according  to  the  interval  {L\.»3)      in  which 
the  corresponding   6.   lies. 


(1^.3) 


A_  J   0  <  ©4  <  1      /  ^-  •  ®i  ^  ^ 


A^  t   -1  <  e,  <  0    I  B^  2  e  <  -1 


-  li^  - 


The  characteristic  slopes,   X.,   and  integers,   j,   may  be 
divided  into  similar  classes  by  the  same  criteria;  we  shall 
use  the  same  notation  for  these  classes. 


B+ 


Pig.  1 


Following  the  above  procedure  we  write   (3.10)   as 


where  for  simplicity  we  have  introduced 


-  15  - 


Since  the  mesh  spacing  is  in  general  "small",  the  inte- 
gral in   (i+.Ij-)   can  be  approximated,  say  by  a  modified 
trapezoidal  rule,  to  yield  the  approximate  equations 

Here   P.   is  some  point  on  the  characteristic  segment  from 
P,   to  Q,   and  from   (3.7),  (4.2)   and  {i\..3) 


We  introduce  quantities  w.  ,   which  are  to  approxi- 
mate  the  corresponding  function  values  (^At    ,r,  )   on  the 
mesh  (i|.0).   At  intermediate  points,  such  as   P.,   we 
approximate  cj>.      by  linear  interpolation: 


(IJ..8)   fi(Pj) 


li^l"l,kil*'l-li^l'<i?l  .   ^/Bj 


There  are  several  obvious  ways  of  approximating  \)r .  ^  /_(P')j 
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for  example  centering  on  the  mesh  diagonal 
(4.9a)   ♦j.i/jCP])  ^\  4  ^  ^ 


or  one  could  center  on  the  characteristics  using   (4.S)  in 


*M/2'^;->  *  ♦j-1/2'^'   • 


Here  Q,^  are  the  points  (t  ,r,  ^,  )  of  Pig.  1  and  vector 
addition  Is  Implied.  The  latter  centering  seems  more  ac- 
curate but  (i4..9a)  Is  simpler  and  we  shall  use  It.  From 
{l\.*5)      we  can  write   (Ij..9a)   more  explicitly  as 

C4.9b)    *j.i/2(i^)  =  -(V  ^)i'"5>"jAn'-<^ic-^)|(w^:i_^. 


Q,  +  4       Q,  +  Q 
"j-l,kil>*Sj(^2-)-^  ^..,(-^2-)  i 


where : 


(i^-.lO) 


^k  -  ^kll/2  ^  i^V  ^k+l^ 
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Introducing  the  approximations   (4.8)-(L|..9a)   into 
(4.6)   and  requiring  the  equality  to  hold  we  obtain:   For 


(""V"j-i,ic'^-D]^*j-i/2'^)  : 


For   P.£B 

J   _   » 

j,lc+l    j-l,k+l^     D.   *j-l/2^   2   ^ 

Inserting   (ij..9b)   Into   (I4..II)   and  collecting  like  terms 
we  may  write : 
For   P. 6  A 

"^•12-'  -T.i  '   ^,,f  l>^""l.K^  "f.-c-l..!^  <=J,."J-1,.11  * 

For  P.6B^ 
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''■''''  ""^ '-  "Lf  J-'^^'J-^''^'  =l^"l^±l^  =],."!-l..ll  * 

J 

The   coefficients   are   defined  by 

j,k  2D.      k     i*]^  J»k        'j      2D.        k     r^n    , 


A4.    _        b, 


^1 


The  source  terms  are  easily  obtained  from  (3.3)   as,  say 


-  ^l/  n+1  ,   n 
1=0 


^^a^j(?^)t(w—  -  w'-^^,^) 


The  Initial  conditions   (3*6)   are  replaced  by 

and  the  boundary  conditions   (3.6)   by 
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(Ij..l6a)  w"?  K  =  °   '    J  <  /2 


Equations  (i|.12)-(I|..l6 )   form  a  determined  system 
for  the  unknowns  w^  ,  .   The  solution  of  this  system  Is 
considered  in  the  next  section.   It  is  clear  that  if  the 
\|f .  ,  /p  had  been  centered  on  the  characteristics  rather 
than  as  in   ([|..9),   equations   ([4.. 12)   would  remain  un- 
changed in  form  while  the  coefficients   c.  ,   of   (I4..13) 
would  be  modified  in  an  obvious  manner.   The  solution  of 
the  resulting  system  could  then  be  analyzed  as  in  the 
next  section. 

The  integral  equations   (3.12)   can  be  approxi- 
mated in  an  analogous  way  yielding  equations  similar  to 
(1|.12),   However  the  coefficients   c.  ,   are  then  approx- 
imations  to  exponential  functions  which  must  be  positive. 
This  procedure,  in  part,  eliminates  difficulties  leading 
to  negative  intensities  which  frequently  arise  when  0" 
or  b  •/„  become  much  larger  than  the  mesh  spacing  (i.e. 
small  mean-free-path  and  origin  difficulties).   More 
care  in  centering  and  approximating  the  integral  term 
on  the  right  of   (3.12)   allows  the  treatment  of  prob- 
lems with  rapid  spatial  and/or  temporal  changes  In 
cross-sections  and  sources.   These  considerations  have 
been  applied  successfully  in  some  classified  problems 
dealing  with  radiation  transport. 
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5«   Solution  of  Difference  Equations  (Converp;ence  of  Iterations) 

In  the  absence  of  scattering  and  fission   (cl   =  0) 
the  system  of  linear  equations   (i|..12)-(i|.l6 )   can  be  so 
ordered  that  the  coefficient  matrix  Is  triangular  and  the 
solution  is  obtained  by  a  simple  recursion.   The  ordering 
and  recursive  solution  Is  simply  described  In  [l].   How- 
ever, If  either  scattering  or  fission  are  present   a. ,f   0, 
and  a  direct  solution  does  not  seem  feasible.   In  these 
cases  an  Iterative  procedure  is  used  which  evaluates   S. 
at  the  old  Iterate  and  solves  for  the  new  Iterate  as  in 
the  non-scattering  case.   We  proceed  to  determine  a  con- 
dition on  the  mesh  spacing  which  is  sufficient  for  the 
convergence  of  this  procedure. 

The  iteration  is  defined  explicitly  by  assuming 

w.  ,   known  for  all   (j,k)  and  some  n  >  0.   Then  we  de- 
i  »^  ~" 

termine  the  iterates,   w.  ,  '^ ,    p  =  0,1,.,,,,   by  taking 
some  initial  guess,  say 


{5,0a)  w^"*"^'^  =  w"" 


and  requiring  each  iterate  to  satisfy  the  bovindary  con- 
ditions  (ij..l6); 


(5.0b)         w^'^J'P  =  0  ,    j  <  V2   , 


'J,K 


(5.0C)  -IIo'^  =  -j:j;g  >  J  >  '/2  . 
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Defining  the  source  iterates  by 


(5.1)  sP-^(-V)  =|[gj(Q.)-gj(^)]  -I-ij^?ic)|(<J'^-'-w?^^,,)  , 


1=0 


we  require  from  (14..  12): 
For   P^€A^ 


(5.2a)  wTi'^=h^     i 


•^■'^    "j,k 


°j,k'^-l',k  ■*■  °j,k''j,k+l  ""  °j,k'^-l,k+l  "" 


'j,k^  j,k   "j-l,k^   2Dj. 


For  P,£  B_|_ 

(^  2b)   W^'^^-^P  =  in    Tc^   w'^'^^'P  +  C^   W^       +  C^   W^         + 

^^•^^^   ^j,k       ^^J^l°j,k''j-l,k  "^  °j,k^J,k+l  ^  °j,k''j-l,k+l 

,  ^ij.   r„n+l,p  ,  „n+l,p   1+  ^iirqP-l/2±!5^+  qP-lfSi^ll  V 

As  previously  indicated  the  system  (5»0)-(5«2)  can  be 
solved  explicitly  for  the  p-th  Iterate  by  the  recursion 
presented  in  [l]. 

The  convergence  of  the  above  procedure  may  be 
examined  by  introducing  the  quantities 

/    »  n+l,p  _  n+l,p    n+l,p-l 
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Then  subtracting  {B»2)      for  the  (p-l)-st   Iterate  from 
that  for  the  p-th  iterate  and  taking  absolute  values  yields 


=1 


(S.^a)   |v]:i'P|<|-J*^J  |v5:j;P|  .  |tP;1|   ,  fo.  P^eA,  , 


c ,  ,     , ,       c^ 


+  |tP;J|,    for   Pj£B^ 


Here  we  have  introduced,  using   (5«1)» 

A  £    J 

Let  us  also  introduce  some  norms  of  the   v.   *^  and  other 
quantities  as  follows 


(5.6a)  f P  =  max  |vf I'Pj   , 

'>      0<k<K  '''^ 


(5.6b)  PP  =  max  f P  =  max  lv^''J:'P| 

■  0<j<J   J   (j,k)   J'^ 

(5.6c)  tP-^  =  max  |tP-J| 

(j,k)   J'^ 


-  23  - 


(5.6d)  p  5  max  -^ 


(5.6e)  a  =  max  |a.  .(?,  )| 


From  the  definitions   (i^..l3)   and   (5.6d)   we  have  since 


(5.7a)  0  <  p  <  1   , 

and 

(S.7b)  max  l-^l  <  1  . 

Taking  maxima  in   (S.i^-)   and  using  the  above  definitions 
and  results  we  obtain 

fP  <  fP_^  +  tP-1  ,   for   jCA^   J 

(5.8) 

fj  <  fj.i  -^  P(fj  +  f^.i)  +  tP-^    ,   for   jGB^   . 

V/ith      (5.7a)      we   can  easily  combine    the   above    into 
(5.9)  fP  <  i±£  rP.,.  Oj  tP-1      ,        for      jeA^DB^      . 

-   2k.  . 


The  sets  A   and  B   have  been  excluded  for  the  moment 
since  at  k  =  0  equations   ($.0c)   rather  than   (5.2) 
hold  for   j^A  UB    Thus  If  the  maxiitium  in  (5.6a)   occurs 
at  k  =  0  we  have,  in  place  of   (5.9), 

f P  <  f P  ,     for   j6  A  UB   and  max  at  k  =  0  . 
J  —  «J-J 

However,  if   j£A  UB   then   j  >  "^ /Z     and  J-jfeA^UB^. 
Thus  by   (5.9) 

(5.10)    fP  <  -^  fP  .+  ^  tP~^  ,    for  j£A  UB  and  max  at  k=0  . 

J  ""  -L-P   «J-J   -L-P  -   -         ^ 

If  this  is  not  the  case  for  some   j   then  the  inequality 
(5.9)   is  valid,  and  we  have  inequalities  for  all   f?. 

J 

Applying   (5«9)   and   (5.10)   recursively  in   J  we  obtain 
for  all   j 


fP 


1  -  P  ^o  ^  i-s  i-p  ^      » 


where 


(5.11)  p  =  fi^ 


>p  _ 


However,  since   (j=0)6;a^UB^  and  f^^=  0  we  have  by   (5.9) 
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f  <  1  tP-1   , 

o  —  1-p        * 


and  thus  from  the  above 


(5.12)  fP  <  ^^ ^  tP"^ 

^^    ''  ^j  -  1-p   1-p 


From  definitions   (5.6b,c,e)   and  {5»$)     we  must 
have 

tP-1<  ^  jaFP-1 

where 

4^5  max  ^   ,  and  by   (i^.7) 

1    A      ^ 

(5.13)  =  max  ^(v^t,  lawpD   , 


V2' 


since   laj/pl  =  ^^^^  |a^|.  Using  the  above  in  (5.12) 

J 

and  taking  the  maximum  with  respect  to   j   yields 


,5.14)       PP  <  [^  |a_  ij^flij  pP-l   . 

Applying   (5«14)   recursively  in  p  implies 

11m  fP  — >  0   , 
p — >oo 

provided  F   is  bounded  and 
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That  this  Implies  convergence  is  easily  seen  from   (5»6b) 
and  (5»3)   by  considering 


tt^  -Id  n+1  _  n+1,0  .    -,  .  f—   „n+l,p 

p'=i 


Condition  (5.l5)   is  thus  sufficient  for  the  con- 
vergence of  the  proposed  iterative  solution.   To  examine 
the  restrictions  it  Imposes  on  the  mesh  spacing  we  note 
from   (5.9)-(5.11)   that  condition   (5.l5)   still  implies 
convergence  if  we  replace   p  and  TZ7,     ^y  their  respec- 
tive bounds,  from   (5»6d)   and   (i|.13)   * 


(5.17a) 


where 


P  <  ^/y  , 
1-p  ^  /^   ' 


a.vAt     Ar 

(5.17b)       Y  5  minYj.  =  min(  \-j^^\ .     I^-^^T^D    • 

Thus,  since   0  <  y  <  1,  (5.l5)  may  be  replaced  by 

A4   ,    p  J+1    -1 

2D  ^  Ja  ^V 

Prom  the   definition   (5.13)   we   obtain,  with  x  =  ^-^/At    , 

Ac 

=  /\t  max  -^^v.  I- 

^j/2 
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^  =  At  i>iax|(v,|^|)      ; 


and  the  above  Inequality  becomes 

-1   o  J+1    -1 

Since  ^An  <  |a.|  <  1  we  have,  by  (5.17b)   , 


(5.18)       At  <  j~t™axi(v,|^|)]      [(|)        -  1] 


(5.19a)        Y  =  min(j4ii   ,^)   = 


'Vv  ,    if     Vv  <  /^V2   , 


and  similarly 

[  V2    ,   if  Vv  <  A^^/2    , 
(5.19b)   inax|(v,|^|)  =< 

VVAt^   ,  if  Vv  >  ^^/2   . 

Using   (5.19)  in  (5.18)  we  obtain,  recalling  that  ^\i  <  2 
and  /Vl  =  /j  ,  three  distinct  forms  for  the  bound  on  At: 

^^^  avJ^2v^  »        V  -     2  * 

(5.20)      (b)     At   <]_2        (^L.)*^""^  4l^<i  <,At 


(c) 


;2   ^2^^  '     ~2"-^  V  -V    2 

V  0_^  »  yj       2.        ~     V 


J+1 


2  2 

Here  we   have   used   the    fact   that     —  >  2     and  hence    (—)        »  1, 

Y  Y 

Conditions  (a)  and  (b)  may  be  applicable  in  practical  com- 
putations while  (c)  must  apply  as  /\v.  — >  0.   In  all  cases 
it  may  be  observed  that,  for  fixed  a,x  and  v.  At  must 
decrease  as   J   increases.   Also,  as  expected,  large 
scattering  and/or  fission  cross-sections  require  a  finer 
mesh  for  convergence.   This  is  easily  seen  from  (5.6e)  and 

(3.3h)  which  yield 

Ja  =  2  max  [vo^(rj^)  +  o^(r^)K(  ti^,H,,)  ]   . 
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6.   Convergence  and  Stability  of  Hyperbolic  Difference  Equations 

In  this  section  we  proceed  to  show  that  the  solution 
of  the  difference  equations  given  by   (i+.H),  (i^..9b)   and 
(i4..1i4--.l6)  converge  as  At  — >  0,  /\r   — >  0  to  the  solu- 
tion of  the  hyperbolic  system   (3.l4-»6),   Independently 
of  the  mesh  ratio,  ^r-  =  x.   Of  course  condition   (5.l5) 
must  be  satisfied  If  scattering  or  fission  Is  present,  but 

-X- 

this  plays  no  part  in  the  present  considerations.   Further 

the  results  of  this  section  are  easily  used  to  show  that 

the  niimerlcal  solution  Is  stable  with  respect  to  small 

changes  of  the  Initial  or  bovmdary  data. 

The  proof  proceeds  by  obtaining  bounds  for  a  system 

of  difference  equations,  closely  related  to  those  of  (i^-.H)^ 

n 


for  the  quantities  v.  ,  : 


For  P  .  £  A^ 


'^•°^'     ^>  '^j'l.^ll-'^-^j'^.k^  ^J*J.l/2'¥^'  *  ^  J.k   • 


For   P.e,B+ 


<^-*'  ^5ik-  ^3^,.ii-'^-^j'^:Jii  =  «j*j-i/2<^'  ^  ^i, 


Here  we  have  defined   y.   as  in   (^.13)  >  and 


(6.0c) 


o,k  -  "^o,k    * 


V^   5  v^   +  v^      : 


Since,  for  fixed  J,  condition  (5.l5)  is  of  course 


guaranteed  in  the  limit  /\t,   /\r   — >  0. 
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the  X  •  and  6.   are  at  present  arbitrary  and  ^.   ,  /p   in 

(6.0)  are  to  be  the  same  runctlons  of  v.  ,  as   ^4  i/p   ^^  ^^ 
the  w^  ,  (as  defined  by  equations  (l4..9b)  and  ([^..11+)).   The 
system  (6.0)  may  be  written  symbolically  for  future  reference  as 

Bovindary  and  initial  equations  similar  to  those  of   (]4..1ij.-.l5) 
are  also  assumed. 

Let  us  introduce  some  norms  by 


(6.2) 


']  =-  T  "^'^^ 

» 

f .  =  max    |v 
^          k 

F^  =  max  F^ 

» 

f  ^  =  max  f  ^ 

For  future  use  we  observe,  from  (6.0c),   the  identity 


'^•2'       ^5,k  =  ^",.  -  ^"-1,.  ^ ^'-I'X,.  • 


Applying   (6.2)   this  yields 

f^  <  ( j+1)  F^   ,  and  thus 
(6.4) 

f^  <  (J+1)  F^    . 


We  now  solve  each  of   (6.0a)   and   (6.0b)   for  V.  ,  ,   take 
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maximuin  absolute  values  with  respect  to  k  and  recall 
that   0  <  Y^  <  1   for   P.€  A^,   0  <  y  <  1   for   P,€  B^ 
to  obtain: 
For   P.^A^ 


(6.sa)        F^^i  <  P^  .  ™ax  |6.*j.,/2(^)l  *  r-    : 


For      P.€  B^ 


Q^^  " 


(6.5b)  p-^l<P-.:nax    M  ,  ._^/^^)  |    .^      . 


Here   we   have    Introduced 

"X  A   =  I'lfLX    It   •   1,1      >      ^"d  for   future   use 
(6.6)  '         "  '' 

rn  /T-n 

=  max   X  . 

J  -^ 

Taking  raaxiina   in      (6.5)      and  applying      (6. 2),  (6. 6)      we   get: 
For      P.e.  A^ 

(6.7a)  F^""^  <  f"  +  max    \6  .^ ^^^^^i-^)  \    +  r  ""      i 


For      P  .€-  B^ 
J       1 


Q^-H^  ^ 


(6.7b)  F^-"^   <  F^   +  max    lIitj.i/2(-|~)  I    +  ^        ; 
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where   y  E  min  y..   ir  the  maxima  In   (6.5)   were  to  occur 
for   k  =  0,   a  consideration  similar  to  that  of  section  5 
easily  shows  that  equations   (6.7)   are  still  valid. 

From  the  definition  of  ^-.n/p  ^^  obtain,  using  {^.•9) 
and  (i|.lij.): 
For   P€:a^UB_^ 

(6.8a)   xnaxU^_^/2(^^l^  ^j^^r^^^j-l''^j"'^j-l^  "^ 


■^      ,n+l^^n, 


i=0 


where 

b. 


aj.ina^dS-J.Iji,,   ,   p,  .  .  ■na.da^jCf,)  |.|a,^  j.,(f,)  |  ) 


(6.eb)  G  =  2  ma>'  tg.-(t  .r,  )| 

j»k,n   *• 


Taking  maxima  with  respect  to   j   in   (6.8a)   gives 


(6.9a)  m&x\^   ^^^^\<   A(f^'^^+  f^)  +  G  <  (  J+1)A(f'^'^^+  F^)  +  G 
where  we  have  used   (6.[|.)   and  the  definition 


( 6 . 9b )  A  =  2  max  a  .  +  max  V"  p .  . 


1=0 
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Using   (6.9a)   in  (6.7)  with 

5 


6  = 
(6.10) 

t""  =  max  (r",  tVr) 


=  max  (16  |,|-i|) 


we  obtain,  provided  A(J+1)5  <  1   , 

(6  n  )       T?^+l  ^  1+A(J+1)6  pH  .   ^"^^6G 
^^•^^^       ^    ^  1-A(J+1)6  ^  "^  1-A(J+1)6 


=n , 


This  implies 

f"  <  B  P  +  Hb"  1-A(J+1)6    » 
(6.12) 

„  _  H-A(J+1)6 

^  =  1-A(J+1)5    * 

and  hence   using      (6.i|.) 


if.   T\)  f"  ^   f  T+1  m^T?°   +  l-B^   (J^H)(f^^6G) 

(6.13)  f      <    (J+1)B   F      +  ilB-  l-A{J+l)6 

This   is   the   desired  bound  on  the      v.   ,       in  terras   of  the 

J  *^ 

initial  valuer   F  ,   given  coefficients  in  A  and  B, 
inhomogeneous  terms  in  ^  ,   and  number  of  equations,  (J+1) 

To  apply  the  above  results  we  note  that  the  dif- 
ference equations   (ij..  11)   may  be  written  as 

(6.1W       l{w^-}=^*^.,,,{.-i] 
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where 


(6.1i^b) 


n 


W"   =  w 


^,k  -  ^5,k  ^  vi,k    • 


If  the  functions   <9.(t,r)   form  a  solution  of  the  hyper- 
bolic system   (3.i4--.6)   with,  say,  second  continuous  de- 
rivatives it  is  easily  shown  that 


(6.15a) 


ID' 


n+1 

J, if 


A5  1  f  n+1^ 

2Df  *j-l/2     jfj,k 


+  J-p  .   , 

2Dj    J,lc 


where 


f j,k  ^     ^j^^n'^k^ 


(6.15b) 


"'0,k   -    ^0,k      * 


^l.  =  <i'5,^  ^  *5-i.i<    ' 


n 


and  P I   -[^      is  the  truncation  error  in  approximating  deriva- 
tives   and  integrals  by  difference  quotients  and  sums.   It 
is  a  tedious  but  elementary  matter  to  show  that 


(6.16) 


Pj^k  =  ^^A*)  +  ^(Ar: 


if  the  coefficients  in   (3.6)   are  bounded  and  have  bounded 
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derivatives.   However  this  is  not  true  due  to  the  — 

r 

term  near  the  origin.   Thus  we  shall  exclude  a  small 
sphere  of  radius  R^  around  this  point,  replace  the  origin 
condition  (3»5b)   by 

^j(t,Rg)  =  0    j  >  V2   , 

sjid  prove  convergence  for  this  altered  problem. 

Let  the  error  between  the  exact  and  numerical  solu- 
tions be 

and  define   v.  ,   as  in  (6.0c).  By  the  linearity  of  the 
operators   L  f*^  and  ^^.i/pf*  \    ^®  have  from  (6.1i|-.l5) 

if  we  set   g.(t,r)  =  0. 

Comparing  this  equation  with   (6.1)   we  see  that  the  bo\md 


(6.13)   applies  if  we  define   G  =  0, 

5  =  max(^^,  P.6A^  ;  TTD'  ^i^^+ 
(6.19a)  j  ^°j    J   -    ^r^i        J   i 


=  V^t 


''The  proof  also  holds  for  a  reflecting  sphere  with  f.(t,R.  )  = 
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and 

Using  the  above  in   (6.13)  we  obtain,  since   (1  +  — )  <  e  , 

2(J+l)AVt    ^    2(J+l)AVt    Ai  Ai-\+M  /\r>) 
(6.20)  |v^  ^|<(J+l).e         ^  F°+[e         n,^^ jglAillfl^i^ 

Thus  for  a  fixed  angular  mesh   |v^   |  ->  0  as   (At.Ar)  ->  0 
provided  the  initial  error,   F°— >  0,   with  the  refinement  of 
the  mesh.   This  implies  convergence  of  the  numerical  solu- 
tion to  the  solution  of  the  hyperbolic  system  of   (J+1) 
equations. 

A  proof  that  the  solution,   tf>.(t,r),   of  the  hyper- 
bolic system  (3.4-»6)  converges  to  the  solution,   (D(t,r,ix), 
of  the  transport  problem  (2.0-,3)  as  J  — >  oo  has  not  been 
completed.   However,  anticipating  this  result,  we  may  assume 


(6.21)       max  |lD(t,r,ti  .)^  d),(  t,r)  |  =  ^(^)   . 
t,r,j        J   ^J  '^ 

Then  using  (6.17)  and  the  above 

Thus  for  convergence  it  is  sufficient  that   |v.  ,  |  =  (9^{-j) , 
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Prom  (6.20)  we  see  that  this  imposes  a  very  stringent 
restriction  on  the  mesh  spacing: 


(6.22)  At  =  &iAr)    <-T^e   2 


where   K,   and  Kp  are  positive  constants.   The  mesh  ratio, 
*-^  //\t ,      is  still  arbitrary.   For  some  steady  state  problems 
the  conditions  on  the  mesh  have  been  shown  to  be  much  less 
restrictive  [2], 

The  above  results  are  easily  used  to  estimate  the 
growth  of  rounding  errors  in  the  nvtmerlcal  solution.   We 
simply  let   v  be  the  error  between  the  exact  and  numerical 
solutions  of  the  difference  equations,  and  replace   F   and 
[  t^(  At )  +  0^{/\v )  ]   in  (6.20)  by  the  maxim\am  rounding  error 
allowed  at  each  stage.   It  is  seen  that  the  computations 
are  stable  with  respect  to  such  errors,  for  fixed  J,   as 
they  are  boumded  independently  of  the  number  of  time  steps, 
n,   and  spatial  steps,   K,   used  to  integrate  up  to  some 

fixed  time,   t  .   So  far  as  one  can  tell  from  our 
n 

bo\and,  an  increase  in  the  number  of  angular  steps,  J, 
might  require  an  exponential  increase  in  the  nvunber  of 
significant  digits  in  order  to  maintain  a  constant 
rounding  error. 
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